Dynamics of an Idealized Model of Microtubule Growth and Catastrophe 



00 
O 

o 

(N 

< 
(N 



0\ 
6 : 
x> : 



(N 
> 

o 
o 
co 
o 
l> 
o 



I , 



X 



T. Antal, 1 P. L. Krapivsky, 2 S. Redner, 2 M. Mailman, 3 and B. Chakraborty 3 

'Program for Evolutionary Dynamics, Harvard University, Cambridge, MA 02138, USA 
Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 02215, USA 
!1 Martin Fisher School of Physics, Brandeis University, Waltham, MA 02454, USA 

We investigate a simple dynamical model of a microtubule that evolves by attachment of guanosine 
triphosphate (GTP) tubulin to its end, irreversible conversion of GTP to guanosine diphosphate 
(GDP) tubulin by hydrolysis, and detachment of GDP at the end of a microtubule. As a function 
of rates of these processes, the microtubule can grow steadily or its length can fluctuate wildly. In 
the regime where detachment can be neglected, we find exact expressions for the tubule and GTP 
cap length distributions, as well as power-law length distributions of GTP and GDP islands. In 
the opposite limit of instantaneous detachment, we find the time between catastrophes, where the 
microtubule shrinks to zero length, and determine the size distribution of avalanches (sequence of 
consecutive GDP detachment events). We obtain the phase diagram for general rates and verify our 
predictions by numerical simulations. 

PACS numbers: 87.16.Ka, 87.17.Aa, 02.50.Ey, 05.40.-a 



I. INTRODUCTION AND MODEL 

Microtubules are polar linear polymers that perform 
major organizational tasks in living cells [3, 0]. Through 
a unique feature of microtubule assembly, termed dy- 
namic instability @ , they function as molecular machines 
[3| that move cellular structures during processes such as 
cell reproduction 0, A surprising feature of micro- 
tubules is that they remain out of equilibrium under fixed 
external conditions and can undergo alternating periods 
of rapid growth and even more rapid shrinking (Fig. [1]). 

These sudden polymerization changes are driven by the 
interplay between several fundamental processes. Micro- 
tubules grow by the attachment of guanosine triphos- 
phate tubulin complexes (GTP) at one end [3|,|6||. Struc- 
tural studies indicate that the end of a microtubule must 
consist of a "cap" of consecutive GTP monomers p} for 
growth to continue jf|. Once polymerized, the GTP of 
this complex can irreversibly hydrolyze into guanosine 
diphosphate (GDP). If all the monomers in the cap con- 
vert to GDP, the microtubule is destabilized and rapid 
shrinkage ensues by the detachment of GDP tubulin 
units. The competition between GTP attachment and 
hydrolysis from GTP to GDP is believed to lead to the 
dynamic instability in which the GTP cap hydrolyzes to 
GDP and then the microtubule rapidly depolymerizes. 
The stochastic attachment of GTP can, however, lead 
to a rescue to the growing phase before the microtubule 
length shrinks to zero [H, Q . 

The origin of this dynamic instability has been ac- 
tively investigated. One avenue of theoretical work on 
this dynamical instability is based on models of mechani- 
cal stability [1, Gil HH • For example, a detailed stochastic 
model of a microtubule that includes all the thirteen con- 
stituent protofilaments has been investigated in Ref. [l(| • 
By using model parameters that were inferred from equi- 
librium statistical physics, VanBuren et al. [lOj found 
some characteristics of microtubule evolution that agreed 
with experimental data [l2j ■ The disadvantage of this de- 



tailed modeling, however, is its complexity, so that it is 
generally not possible to develop an intuitive understand- 
ing of microtubule evolution. 

Another approach for modeling the dynamics of mi- 
crotubules is based on effective two-state models that 
describe the dynamics in terms of a switching between 
a growing and a shrinking state d, [HI, 0, [H, [H, [TtJ • 
The essence of many of these models is that a micro- 
tubule exists either in a growing phase (where a GTP 
cap exists at the end of the microtubule) or a shrinking 
phase (without a GTP cap) , and that there are stochastic 
transitions between these two states. By tuning param- 
eters appropriately, it is possible to reproduce the phase 
changes between the growing and shrinking phases of mi- 
crotubules that have been observed experimentally p|. 
While the two-state model has the advantage of having 
only a few parameters, a constant rate of switching be- 
tween a growing and shrinking microtubule is built into 
the model. Thus switching models cannot account for 
the stochastic avalanches and catastrophes that occur in 
real microtubules. 

On the other hand, a minimalist model of microtubule 
dynamics has been proposed and investigated by Fly- 
vbjerg et al. In their model, they dispense with 

attempts to capture all of the myriad of experimental 
parameters within a detailed model, but instead con- 
structed an effective continuous theory to describe mi- 
crotubule dynamics. Their goal was to construct an ef- 
fective theory that contained as few details as possible. 
As stated in Ref. [Hj], they envision that their effective 
theory should be derivable from a fundamental, micro- 
scopic theory and its parameters. 

This minimalist modeling is the approach that we 
adopt in t he p resent work. We investigate a recently 
introduced [TiJ [2(| kinetic model that accounts for many 
aspects of microtubule evolution. Our main result is that 
only a few essential parameters with simple physical in- 
terpretations are needed to describe the rich features of 
microtubule growth, catastrophes, and rescues [2lj . 
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We treat a microtubule as a linear polymer that con- 
sists of GTP or GDP monomers that we denote as + 
and — , respectively. To emphasize this connection be- 
tween chemistry and the model, we will write the former 
as GTP + and the latter as GDP - . The state of a micro- 
tubule evolves due to the following three processes: 

1. Attachment: A microtubule grows by attachment 
of a guanosine triphosphate (GTP+) monomer. 

| h) =H h +) rate A 

| } =H h) rate pX. 

2. Conversion: Once part of the microtubule, each 
GTP+ can independently convert by hydrolysis to 
a guanosine diphosphate (GDP - ). 



+ 



rate 1. 



3. Detachment: a microtubule shrinks due to detach- 
ment of a GDP - monomer only from the end of the 
microtubule. 



>=>!■■■) 



rate \x. 



Here the symbols | and ) denote the terminal and the ac- 
tive end of the microtubule. It is worth mentioning that 
these steps are similar to those in a recently-introduced 
model of DNA sequence evolution [22j , and that some of 
the results about the structure of DNA sequences seem 
to be related to our results about island size distributions 
in microtubules. 

Generically, the (X,fi,p) phase space separates into a 
region where the microtubule grows (on average) with a 
certain rate V(X,/j,,p), and a compact phase where the 
average microtubule length is finite. These two phases 
are separated by a phase boundary fx = /i»(A,p) along 
which the growth rate V(A, fj,,p) vanishes. While the be- 
havior of a microtubule for general parameter values is 
of interest, we will primarily focus on extreme values of 
the governing parameters where we can obtain a detailed 
statistical characterization of the microtubule structure. 
For certain properties, such as the shape of the phase 
diagram, we will also present results of numerical simu- 
lations of the model. 

In Sec. UH we study the evolution of a microtubule un- 
der unrestricted growth conditions — namely no detach- 
ment and an attachment rate that does not depend on 
the identity of the last monomer. Our results here are 
relevant to understanding the distribution of cap length 
and the diffusion coefficient of the tip of the microtubule 
in the growth phase. The predictions of the model in this 
limit could also be useful in understanding the binding 
pattern of proteins to microtubules [23| . Since proteins 
are important regulatory factors in microtubule polymer- 
ization, these results could prove useful in interpreting 
the effects of proteins on microtubule growth. 

By a master equation approach, we will determine both 
the number of GTP + monomers on a microtubule, as well 
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FIG. 1: Numerical simulations of typical microtubule lengths 
versus time for detachment rate fi — 5 and attachment rates: 

(a) A = 1.4, where the microtubule generally remains short, 

(b) A = 1.5, where the length fluctuates strongly, and (c) 
A = 1.6, where the microtubule grows nearly steadily. 



as the length distributions of GTP + and GDP - islands 
(Fig. Many of these analytical predictions are veri- 
fied by numerical simulations. In Sec. IIII[ we extend our 
approach to the case of constrained growth, p / 1, in 
which microtubule growth depends on whether the last 
monomer is a GTP + or a GDP - . In Sec. IIV1 we inves- 
tigate the phenomenon of "catastrophe" for infinite de- 
tachment rate /i, in which a microtubule shrinks to zero 
length when all of its constituent monomers convert to 
GDP - . We derive the asymptotic behavior of the catas- 
trophe probability by expressing it as an infinite prod- 
uct and recognizing the connection of this product with 
modular functions. We also determine the asymptotic 
behavior of the size distribution of avalanches, namely, 
sequences of consecutive GDP - detachment events. Fi- 
nally, in Sec. El we discuss the behavior of a microtubule 
for general parameter values through a combination of 
numerical and analytic results. Here numerical simula- 
tions are useful to extract quantitative results for param- 
eter values that are note amenable to theoretical analysis. 
Several calculational details are given in the appendixes. 
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II. UNRESTRICTED GROWTH 

We define unrestricted growth as the limit of detach- 
ment rate /z = 0, so that a microtubule grows without 
bound. Here we consider the special case where the at- 
tachment rate does not depend on the identity of the last 
monomer; that is, the limit of p = 1, where the attach- 
ment is unconstrained. Because of the latter condition, 
the number N of GTP + monomers decouples from the 
number of GDP - , a greatly simplifying feature. 



<L>=A1 



FIG. 2: (Color online) Cartoon of a microtubule in unre- 
stricted growth. Regions of GTP + are shown dark (blue) and 
regions of GDP" are light (yellow). The GTP + regions get 
shorter further from the tip that advances as Xt, while the 
GDP - regions get longer. 



whose solution is an arbitrary function of t — y or, equiva- 
lently, (1 — z)e~ t . If the system initially is a microtubule 
of zero length, IIjv(i = 0) = <5jv,o, the initial generating 
function n(z,i = 0) = 1, so that Q = e -Az = e A ( 1- ^e -A . 
Thus for t > 0, Q = e A(1 " 2)e " t e" A , from which 



U(z,t) = e 



_ -A(l- Z )(l-e -t ) 



(0) 



Expanding this expression in a power series in z, the 
probability for the system to contain N GTP + monomers 
is the time-dependent Poisson distribution 



Iljv(t) = 



_ [A(l-e- t r B -A(i-0 



(7) 



From this result, the mean number of GTP + monomers 
and its variance are 



(N) = (N 2 ) - (N) 2 = A(l-e~*). 



(8) 



B. Tubule Length Distributions 



Distribution of Positive Monomers 



The average number of GTP + monomers evolves as 



Jt (N)=X-(N). 



(1) 



The gain term accounts for the adsorption of a GTP+ 
at rate A, while the loss term accounts for the conver- 
sion events GTP+ — > GDP - , each of which occurs with 
rate 1. Thus (N) approaches its stationary value of A 
exponentially quickly, 



(N)=\(l-e- t ). 



(2) 



More generally, consider the probability II/v(t) that 
there are N GTP+ monomers at time t. This probability 
evolves according to 



dn 



N 



dt 



-(N + X)U N +XIL N _ 1 + (N+l)U N+1 . (3) 



The loss term (N + \)Hn accounts for conversion events 
GTP+ -> GDP - that occur with total rate N, and the 
attachment of a GTP + at the end of the microtubule of 
length N with rate A. The gain terms can be explained 
similarly. 



In terms of generating function Tl(z) = X)w=o ^ 
Eq. ((3]) can be recast as the differential equation 



N 



dli 
~dt 



(1-z) 



dn 

dz 



An 



(4) 



Introducing Q = lie Az and y = log(l — z), we transform 
Eq. ([!]) into the wave equation 



dQ dQ 
dt dy 



0. 



(5) 



The length distribution P(L,t) of the microtubule 
evolves according to the master equation 



dP(L,t) 
dt 



For the initial condition P(L, 0) 
again the Poisson distribution 



= X[P(L-l,t)-P(L,t)} (9) 
Sl.Oi t ne solution is 



P(L,t) 



(Xt) 1 
LI 



-XI 



from which the average and the variance are 
(L)=Xt, (L 2 ) - (L) 2 = Xt . 



(10) 



(11) 



Thus the growth rate of the microtubule and the diffusion 
coefficient of the tip are 



V = X, D = X/2 



(12) 



A more comprehensive description is provided by 
the joint distribution P(L, N, t) that a microtubule has 
length L and contains N GTP+ monomers at time t. 
This distribution evolves as 



dP{L,N) 
di 



= XP(L - 1, N - 1) - (N + X)P(L, N) 
+ (N + l)P(L,N+l). (13) 



This joint distribution does not factorize, that is, 
P(L,N,t) ^ P(L,t)U N (t), because (LN) ^ (L)(N). To 
demonstrate this inequality, we compute (LN) by multi- 
plying Eq. (fT3"|) by LN and summing over all L > N > 
to give 



dt 



(LN) = X((L) + (N) + 1) - (LN). 



(14) 
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FIG. 3: Representative configuration of a microtubule, with a 
GTP+ cap of length 4, then three GTP+ islands of lengths 1, 
3, and 2, and three GDP~ islands of lengths 3, 2, and a "tail" 
of length 5. The rest of the microtubule consists of GDP". 



Using Eqs. 
we obtain 



and (fTTj) for (N) and (L) and integrating 



(LN) = A 2 i(l-< 
= (L)(N)- 



- t ) + X(l-e- t ) 
(N). 



(15) 



Using Eq. (JUJ), we have (LN) = (L)(N)(1 + so that 
the joint distribution is factorizable asymptotically. For 
completeness, we give the full solution for P(L, N, t) in 
Appendix [A] 



C. Cap Length Distribution 

Because of the conversion process GTP+ — > GDP", 
the tip of the microtubule is comprised predominantly 
of GTP + , while the tail exclusively consists of GDP". 
The region from the tip until the first GDP" monomer is 
known as the cap (Fig. [3]) and it plays a fundamental role 
in microtubule function. We now use the master equation 
approach to determine the cap length distribution. 




FIG. 4: Cap length distribution obtained from simulations 
at fi = 0, A = 100, and p = 1 compared to the theoretical 
prediction of Eq. P0) . 

Consider a cap of length k. Its length increases by 1 
due to the attachment of a GTP+ at rate A. The con- 
version of any GTP+ into a GDP" at rate 1 reduces the 
cap length from k to an arbitrary value s < k. These 
processes lead to the following master equation for the 
probability nk that the cap length equals k: 



fik = A(n fe _i - n k ) - kn k + 



s>k+l 



(16) 



Equation ([16]) is also valid for k = if we set n_i = 0. 
Note that n = Prob{— )} is the probability for a cap of 
length zero. We now solve for the stationary distribution 
by summing the first k — 1 of Eqs. (|16p with fik set to 
zero to obtain 



rik-i 



s>fc 



The cumulative distribution, Nk 
lies the recursion 



>fe 



(17) 



thus satis- 



Nh 



A 



k + A 



JV fc _i. 



(18) 



Using the normalization Ao = 1 and iterating, we obtain 
the solution in terms of the Gamma function 1241: 



N k = 



x k r(i + A) 
r(fc + i + A) 

Hence the cap length distribution is 

r(i+A) 
rik ~ r(fc+2+A) 

and the first few terms are 
1 



(ft + l)A* 



(19) 



(20) 



n 



l + A 



2A 



n 2 = 



(1 + A)(2 + A) 
3A 2 

(1 + A)(2 + A)(3 + A)' 



Results of direct simulation of the kinetic model are com- 
pared to the predicted cap length distribution (Eq. ([201) ) 
in Fig.[U Because of the finite length of the simulated mi- 
crotubule, there is a largest cap length that is accessible 
numerically. Aside from this limitation, the simulations 
results are in agreement with theoretical predictions. 

It is instructive to determine the dependence of the 
average cap length (k) — J2k>o ^ Uk on ^ Using nk = 
Nk — Nk+i, we rearrange (k) into 



(22) 



k>l 



Using (j!9p . the above sum may be written in terms of 
the confluent hypergeometric series (24|: 



(k) 



-1 + ^(1; l + A; A). 



(23) 



We now determine the asymptotic behavior of (k) by 
using the integral representation 



F(a; b; z) — 



T(b) 



T(b-a)T(a) J 
to recast the average cap length ([2"3"[) as 

A 



dte zt t^il-t) 



fc-a-l 



(k) 



-1 



+ A(£) 7 (A,A), 



(24) 
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where 7(0, x) — f$ dtt a ~ 1 e~ t is the (lower) incomplete 
gamma function. 

In the realistic limit of A ^> 1, we use the large A 
asymptotics 



to give 



(fc) -> v/ttA/2 as A 



00 



(25) 



Thus even though the number of GTP + monomers equals 
A, only \/A of them comprise the microtubule cap, as 
qualitatively illustrated in Fig. [2l Note that the average 
cap length is proportional to the square-root of the ve- 
locity; essentially the same result was obtained from the 
coarse-grained theory of Flyvbjerg et al. fl8j |. 



D. Island Size Distributions 

At a finer level of resolution, we determine the distri- 
bution of island sizes at the tip of a microtubule (Fig. [3]). 
A simple characteristic of this population is the aver- 
age number I of GTP+ islands. If all GTP+ islands 
were approximately as long as the cap, we would have 
/ ~ (N)/(k) ~ VJ. As we shall see, however, I scales 
linearly with A because most islands are short. A similar 
dichotomy arises for negative islands. 

To write the master equation for the average number of 
islands, note that the conversion GTP + — > GDP~ elim- 
inates islands of size 1. Additionally, an island of size 
k > 3 splits into two daughter islands, and hence the 
number of islands increases by one, if conversion occurs 
at any one of the fc — 2 in the interior of an island as 
illustrated below: 



fc-2 



Conversely, if the cap has length 0, attachment creates 
a new cap of length 1 at rate A. The net result of these 
processes is encoded in the rate equation 



dl 



^(fc-2)/ fe + An , 



(26) 



fe>i 



with Ik the average number of GTP+ islands of size fc. 

We now use the sum rules I — X)fc>i^fe and (N) z 
J2k>i ^ k to recas t ((23) as 



dl 
~dl 



= (N) -21 + Xn 



(27) 



from which the steady-state average number of islands is 



I = 



(N) + Xn _ A 2 + A 
2 ~ 2 1 + A' 



(28) 



For large A, the number of islands approaches A/2, while 
the number of GTP+ monomers equals A. Thus the 
typical island size is 2. Nevertheless, as we now show, 
the GTP + and GDP~ island distributions actually have 
power-law tails, with different exponents for each species. 

The GTP + island size distribution evolves according 
to the master equation 



I k = -kh 



2 E I' 

s>k+l 



\(nk- 



Tlfc) 



(29) 



This equation is similar in spirit to Eq. (I16|) for the cap 
length distribution. As a useful self-consistency check, 
the sum of Eqs. {22]) gives ([27]l . while multiplying ([29]) 
by k and summing over all fc > 1 gives Eq. (fT]). 




FIG. 5: Simulation results at /1 = 0, A = 100, and p = 1 for 
the size distribution of positive islands, Jfe/A. The solid line 
is the theoretical prediction of Eq. (|33p . 



(30) 



The stationary distribution satisfies 

fc/fe = 2 E !s + H n k-i - n k ). 



s>k+l 



Using J2 S >2 Is = I — h, w transform ([30]) at fc = 1 to 

3/i = 21 + A(n - rn) 

Similarly, using X) s >3 ^ = I — h — I2 we transform (|3CT 
at fc = 2 to 

47 2 = 2(1 - h) + A(m - n 2 ) 
Thus using (|20|) and (|30[) we obtain 



h 



A 



A 



4A 



3(1 + A)(2 + A) 



h = — 



25A 2 - 6A 



12 12(1 + A)(2 + A)(3 + A) 



(31a) 
(31b) 



The same procedure gives Ik for larger fc. 

Since the Ik represent the average number of islands of 
size fc, they become meaningful only for A — > 00 where an 
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appreciable number of such islands exist. In this limit, 
we write Ik more compactly by first rearranging (|30|) into 
the equivalent form 

(fc - l)J fc _i - (fc + 2)4 = A(n fc _ 2 - 2n fe _x + n fc ). (32) 



We then use (|20|) and the asymptotic properties of the 
Gamma function to find that the right-hand side of 
Eq. ([32D is 



A(rtfc-2 - 2n fc _i + n fe J = - h C I 



and is therefore negligible in the large-A limit. Thus (|32|) 
reduces to (fc — l)/fc-i = (fc + 2)/fc, with solution Ik = 
A/[k(k+l)(k+2)]. We find the amplitude A by matching 
with the exact result, Eq. (|31a[) , to give I\ — A/3 for large 
A. The final result is 



2A 



k(k+ 1)0 + 2) 



(33) 



In the large A limit, / = A/ 2, and the above result can 
be re-written as 



h 
I 



k(k + l)(k + 2) 



(34) 



Remarkably, the size distribution of the positive islands is 
identical to the degree distribution of a growing network 
with strictly linear preferential attachment [25TI261. |27| . 

The results for the island size distribution in the large A 
limit are compared to simulation results in Fig. These 
asymptotic results are expected to apply to island sizes 
k much smaller than the size of the cap which scales as 
VA. The distributions obtained from the numerical sim- 
ulations should then obey the theoretical form but with a 
finite-size cutoff. The results in Fig. [5] are consistent with 
this picture but, interestingly, the numerical distribution 
rises above the theoretical curve before falling sharply 
below it. This anomaly occurs in many heterogeneous 
growing network models, and it can be fully character- 
ized in terms of finite-size effects [23 . 



E. Continuum Limit, A — > 00 

When A — > 00, both the length of the cap and the 
length of the region that contains GTP+ become large. 
In this limit, the results from the discrete master equation 
can be expressed much more elegantly and completely by 
a continuum approach. The fundamental feature is that 
the conversion process GTP + — > GDP" occurs indepen- 
dently for each monomer. Since the residence time of 
each monomer increases linearly with distance from the 
tip, the probability that a GTP + does not convert decays 
exponentially with distance from the tip. This fact alone 
is sufficient to derive all the island distributions. 

Consider first the length I of the populated region 
(Fig. [3]). For a GTP + that is a distance x from the tip, 



its residence time is r = x/X in the limit of large A. 
Thus the probability that this GTP+ does not convert 
is e~ T = e _a A. We thus estimate I from the extremal 
criterion [29f 



(35) 



X>1 



that merely states that there is of the order of a single 
GTP + further than a distance I from the tip. Since (1 — 
e~ 1//A ) -1 — » A when A is large, the length of the active 
region scales as 

i = A In A (36) 
The probability that the cap has length k is given by 



(l-e-C*+i)A)jj, 



The product ensures that all monomers between the tip 
and a distance k from the tip are GTP+, while the pref- 
actor gives the probability that a monomer is a distance 
k + 1 from the tip is a GDP - . Expanding the prefactor 
for large A and rewriting the product as the sum in the 
exponent, we obtain 



k + l 



+ 1 „-fc(fc+X)/2A 

A 



(37) 



a result that also can be obtained by taking the large-A 
limit of the exact result for given in Eq. ([2H)l . 

Similarly, the probability to find a positive island of 
length k that occupies sites x + l,x + 2, ...,x + k is 

k 

(1 - e- x ' x ){l - e -(*+*+i)A) Y[ e- {x+l)/x . (38) 

j=i 

The two prefactors ensure that sites x and x + k + 1 
consist of GDP", while the product ensures that all sites 
between x + 1 and x + k are GTP + . 

Most islands are far from the tip and they are relatively 
short, k ^ x, so that (13"5)) simplifies to 



(39) 



(1 - e - x/x ) 2 e~ kx/x e" fe2/2A . 



The total number of islands of length k is obtained by 
summing the island density l|39p over all x. Since A > 1, 
we replace the summation by integration and obtain 



h 



dx{\-e- x ' x fe- kx ' x e- k2 ' 2X 



2A 



k(k+l)(k+2) 



-k 2 /2\ 



(40) 



The power law tail agrees with Eq. (|33p , whose derivation 
explicitly invoked the A — > 00 limit. 
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We can also obtain the density of negative islands in 
this continuum description, a result that seems impossi- 
ble to derive by a microscopic master equation descrip- 
tion. In parallel with (|39[) . the density of negative islands 
of length k <C x with one end at x is given by 



-*»A(i_ e -« 



/X\k 



(41) 



and the total number of negative islands of length k is 

Wo dxe ^ ,X{l ~ e ~ x,X)k= (k + i)ik + 2y ^ 

Again, we find a power-law tail for the GDP - island size 
distribution, but with exponent 2. The total number 
of GDP - monomers within the populated zone is then 
Sfe>i kJk- While this sum formally diverges, we use the 
upper size cutoff, k* ~ A to obtain J2k>i ^Jk — A In A. 
Since the length of the populated zone I ~ A In A, this 
zone therefore predominantly consists of GDP - islands. 

In analogy with the cap, consider now the "tail" — 
the last island of GDP - within the populated zone (see 
Fig. [3]) . The probability nik that it has length k is 



-£/\\k 



(43) 



Using 



we simplify the above expression to 
m k = A -1 (l-A -1 ) fc = A -1 e -fe /\ 



Hence the average length of the tail is 



(k) = ^ km k = A, 



(44) 



fe>i 



which is much longer (on average) than the cap. 



III. CONSTRAINED GROWTH 

When p 7^ 1, the rate of attachment depends on the 
state of the tip of the microtubule — attachment to a 
GTP + occurs with rate A while attachment to a GDP - 
occurs with rate pX. While this state dependence makes 
the master equation description for the properties of the 
tubule more complicated, qualitative features about the 
structure of the populated zone are the same as those in 
the case p = 1. In this section, we outline some of the 
basic features of the populated zone when p ^ 1 , but we 
still keep u — 0. 



A. Distribution of GTP+ 

The average number of GTP + monomers now evolves 
according to the rate equation 



dt 



(N) 



-(N) +pXn + X(l-n ), 



(45) 



which should be compared to the rate equation Eq. (TTJ 
for the case p = 1. The loss term on the right-hand 
side describes the conversion GTP + ^GDP~, while the 
remaining terms represent gain due to attachment to a 
GTP + with rate A and to a GDP - with rate pX. Here 
no is the probability for a cap of length zero, that is, the 
last site is a GDP - . The stationary solution to (|45|) is 



(N) = pXn + X{l-n ), 



(46) 



so we need to determine tiq. By extending Eq. (fT6|) to 
the case p ^ 1, we then find that no is governed by the 
rate equation 



ho = -pXn + (1 - n ). 



(47) 



Thus asymptotically no 



1+pA 



and substituting into 



46]) . the average number of GTP+ monomers is 



(AO = pX 



1 + A 
1+pA 



(48) 



More generally, we can determine the distribution of 
the number of GTP+ monomers; the details of this cal- 
culation are presented in Appendix iBl 



B. Growth Rate and Diffusion Coefficient 

The growth rate of a microtubule equals pX when the 
cap length is zero and to A otherwise. Therefore 



V(p, A) = pXn + A(l - n ) = pX 



1 + A 
1+pA 



(49) 



For the diffusion coefficient of the tip of a microtubule, we 
need its mean-square position. As in the case p = 1, it is 
convenient to determine the probability distribution for 
the tip position. Thus we introduce X(L, t) and Y(L,t), 
the probabilities that the microtubule length equals L 
and the last monomer is a GTP+ or a GDP - , respec- 
tively. These probabilities satisfy 



dX{L) 

dt 
dY{L) 

dt 



XX{L-1) +pXY(L-l) - (1 + X)X(L) (50a) 
X(L) — pXY(L), (50b) 



Summing these equations, the length distribution 
P(L) = X(L) + Y(L) satisfies 



dP(L) 
dt 



XX (L - 1) +pXY{L -l)-XX(L)- pXY(L). 

(51) 

The state of the last monomer does not depend on the 
microtubule length L for large L. Thus asymptotically 



X(L) = (1 - n )P(L), Y(L) = n P(L). 



(52) 



Substituting (|52p into (|51[) we obtain a master equation 
for the tubule length distribution of the same form as 
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Eq. ([9]), but with prefactor V given by (|49|) instead of 
A. As a result of this correspondence, we infer that the 
diffusion coefficient is one-half of the growth rate, 



(53) 



For large A both the growth rate of the tip of the mi- 
crotubule and its diffusion coefficient approach the cor- 
responding expressions in Eq. (fT2f for the case p = 1. 



C. Cap Length Distribution 

The master equations for the cap length distribution 
are the same as in the p = 1 case when k > 2. The 
master equations for k — and k — 1 are slightly changed 
to account for the different rates at which attachment 
occurs at a GDP~ monomer: 

pXno = Ni = 1 — no 
(1 + A)«i - p\n — N 2 — 1 — n - m 

Solving iteratively we recover no = i+ p \ and also obtain 

2pX 



m = 



n 2 



(l+pA)(2 + A) 
3pA 2 

(l+pA)(2 + A)(3 + A)' 



(54a) 
(54b) 



etc. The general solution for the n^ is found by the same 
method as in Sec. Iff CI to be 



n k = (k + 1)A* 



P 



T(2 + A) 



1 + pX T(k + 2 + A) ' 



(55) 



which are valid for k > 1. With this length distribution, 
the average cap length is then 



(56) 



and the sum can again be expressed in terms of hyperge- 
ometric series as in Eq. (f2"3"]) . Rather than following this 
path, we focus on the most interesting limit of large A. 
Then the cap length distribution (|55[) approaches to the 
previous solution ([2H)l for the case p = 1 and the mean 
length reduces to ([25]) . independent of p. 



D. Island Size Distribution 

For the distribution of island sizes, the master equation 
remains the same as in the p = 1 case when k > 2. 
However, when k — 1 , the master equation becomes 



I\ = 2 I s + pXno — Xn\ 



(57) 



in the stationary state. Then the average number of is- 
lands and the average number of islands of size 1 are 
found from 



21 
3/i 



(jV) + pXn 

21 + pXno — Xrii 



Using no 



jq^j and Eqs. |48|) and (|54a[) we obtain 
pX 2 + A 



I = 



h 



2 1+pX 



pX 



pX 



A/3 



2 + A 



(58a) 
(58b) 



Again, in the limit of large A, the average island size 
distribution reduces to our previously-quoted results in 
(f33|) or equivalently (|34[) . The leading behavior in the 
A — > oo limit is again independent of p. 



IV. INSTANTANEOUS DETACHMENT 

For /i > 0, a microtubule can recede if its tip con- 
sists of GDP - . The competition between this recession 
and growth by the attachment of GTP + leads to a rich 
dynamics in which the microtubule length can fluctuate 
wildly under steady conditions. In this section, we focus 
on the limiting case of infinite detachment rate, \i = oo. 
In this limit, any GDP~ monomcr(s) at the tip of a mi- 
crotubule are immediately removed. Thus the the tip is 
always a GTP+; this means that the parameter p be- 
come immaterial. Finally, for fi = oo, we also require 
the growth rate A — > oo to have a microtubule with an 
appreciable length. This is the limit considered below. 

As soon as the last monomer of the tubule changes 
from a GTP + to a GDP~ , a string of k contiguous GDP~ 
monomers exist at the tip and they detach immediately. 
We term such an event an avalanche of size k. We now 
investigate the statistical properties of these avalanches. 



A. Catastrophes 

The switches from a growing to a shrinking state of 
a microtubule are called catastrophes [8J]. If a newly- 
attached GTP + at the tip converts to a GDP~ and 
the rest of the microtubule consists only of GDP~ at 
that moment, the microtubule instantaneously shrinks to 
zero length, a phenomenon that can be termed a global 
catastrophe. We now determine the probability for such 
a catastrophe to occur. Formally, the probability of a 
global catastrophe is 



e(A) = 



i 



l + A 



IK 1 



,-n/\\ 



(59) 



s>2 



The factor (I + A) 1 gives the probability that the 
monomer at the tip converts to a GDP~ before the next 
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attachment event, while the product gives the probabil- 
ity that the rest of the microtubule consists of GDP - . 
In principle, the upper limit in the product is set by the 
microtubule length. However, for n > A, each factor in 
the product is close to 1 and the error made in extending 
the product to infinity is small. The expression within 
the product is obtained under the assumption that the 
tubule grows steadily between these complete catastro- 
phes and the smaller, local catastrophes, are therefore 
ignored in this calculation. 

The leading asymptotic behavior of the infinite prod- 
uct in (|59[) is found by expressing it in terms of the 
Dedekind r\ function [301 ] 

00 

r](z) = e l7Tz/12 e 27Tinz ), (60) 

n=l 

and using a remarkable identity satisfied by this function, 

?y(— 1/z) = v '—izr)(z). 

For our purposes, we define a — —iirz to rewrite this 
identity as 



JJ (1 _ e -2a„ )= /ZL e («-*>)/12-Q (1 



-2bn\ 



where b 
yields 



7r 2 /a. Specializing to the case a = (2A)~ 



e(A) = 



1 + A 

f 27T 



-it 2 \/6 „1/24A 



II a 



n>l 



-7T 2 A/6 



(61) 



Since the time between catastrophes scales as the inverse 
of the occurrence probability, this inter-event time be- 
comes very long for large A. 



B. Avalanche Size Distribution 

In the instantaneous detachment limit, ji = 00, the 
catastrophes are avalanches whose size is determined by 
the number of GDP~s between the tip and the first 
GTP + island. A global catastrophe is an avalanche of 
size equal to the length of the tubule, whose occurrence 
probability was calculated in the preceding section. Sim- 
ilar arguments can be used to calculate the size distribu- 
tion of the smaller avalanches. 

Since the cap is large when A is large, an avalanche of 
size 1 arises only through the reaction scheme 



where the first step occurs at rate 1 and the second step 
is instantaneous. Since attachment proceeds with rate 



A, the probability that conversion occurs before attach- 
ment is A\ = (1 + A) -1 ; this expression gives the relative 
frequency of avalanches of sizes > 1. Analogously, an 
avalanche of size two is formed by the events 



-+> 



and the probability that the first two steps occur be- 
fore an attachment event is A2 ~ A~ 2 to lowest order. 
At this level of approximation, the relative frequency of 
avalanches of size equal to 1 is Ai — A2 ~ A -1 . Since we 
are interested in the regime A > 1 , we shall only consider 
the leading term in the avalanche size distribution. 

Generally an avalanche of size k is formed if the system 
starts in the configuration 



fc-i 



then k — 1 contiguous GTP+ monomers next to the tip 
convert, and finally the GTP + at the tip converts to 
GDP - before the next attachment event. The probabil- 
ity for the first k — 1 conversion events is A"'*"- 1 ) (k — 1)!, 
where the factorial arises because the order of these steps 
is irrelevant. The probability of the last step is A -1 . Thus 
the relative frequency of avalanches of size k is 



A k ~ \- k T(k). 



(62) 



The result can also be derived by the approach of 
Sec. Ill El We use the fact that the configuration 



fc-i 



occurs with probability Hi<n<fc 

_ x (l - e- n ' x ). Multi- 
plying by the probability that the monomer at the tip 
converts before the next attachment event then gives the 
probability for an avalanche of size k: 



k-l 



A k = (l + \)- 1 Y[(l-e~ n / x ) 



(63) 



Using 1 — e~ n / x = n/Xwe recover (|6"2")l . If we expand the 
exponent to the next order, 1 — e~"/ A « n/A — n 2 / (2A 2 ), 
Eq. (|63|) becomes 



A k = x- k T(k)Y[ UA- fc r(fc) e 



-fc 2 /4A 



V. GENERAL GROWTH CONDITIONS 

The general situation where the attachment and de- 
tachment rates, A and /1, have arbitrary values, and where 
the parameter p / 1 seems analytically intractable be- 
cause the master equations for basic observables are cou- 
pled to an infinite hierarchy of equations to higher-order 
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variables. For example, the quantity uq = Prob{— )}, 
the probability for a cap of length zero, satisfies the ex- 
act equation 



n 



-pXn Q + (1 - n ) - mNq 



and the speed of the tip is 



V(A, v,p) — pXn + A(l - n ) — /zn 



(64) 



(65) 



Here 3\T = Prob{H — )} is the probability that there is 
a GDP~ at the front position with a GTP+ on its left. 
Thus to determine uq we must find Nq, which then re- 
quires higher-order correlation functions, etc. This hier- 
archical nature prevents an exact analysis and we turn 
to approximate approaches to map out the behavior in 
different regions of the parameter space. 



A. Limiting Cases 

For A 7( u < 1, the conversion GTP+ -> GDP" at 
rate one greatly exceeds the rates X,pX,/j, of the other 
three basic processes that govern microtubule dynam- 
ics. Hence we can assume that conversion is instanta- 
neous. Consequently, the end of a microtubule consists 
of a string of GDP~, | • • • ), in which the tip ad- 
vances at rate pX and retreats at rate /i. Thus from (|65|) 
the speed of the tip is 



V(p, X, n) = pX - [i 



(66) 



when pX> ji. 

On the other hand, for A > 1, no = Prob{— )} is 

small and Prob{ )} is exceedingly small. Hence no = 

Prob{ )}+N ~ No- Substituting this result into (JM]) 

and solving for n we find 



1 



n 



1+pX + fi 



(67) 



Note that indeed n <C 1 when A > 1. Using ([57]) in (ESI) 
we obtain the general result for the growth velocity 



V = X 



(1-P)A + 



1 + pX + fi 



when A 3> 1. 



(68) 



B. The Phase Boundary 

A basic characteristic of microtubule dynamics is the 
phase boundary in the parameter space that separates 
the region where the microtubule grows without bound 
and a region where the mean microtubule length remains 
finite. For small A, this boundary is found from setting 
V = in Eq. (f66|) to give /x» = pX for A« 1. The phase 
boundary is a straight line in this limit, but for larger A 
the boundary is a convex function of A (see Fig. [5]). We 
can compute the velocity to second order in \i and A by 



assuming No = and then using (|64|) and (|65|) . This 
leads to the phase boundary 



^* = pX + pX , 



(69) 



that is both convex and more precise. On this phase 
boundary, the average tubule length grows as y/i. 




FIG. 6: Phase diagrams of a microtubule from simulations 
for (a) p = 1 and (b) p = 0.1. The dashed line represents the 
prediction (|69|) that is appropriate for small fj,. The extremes 
of the error bars are the points for which the tubule velocities 
are 0.005 and 0.015, and their average defines the data point. 



When A is large, Eq. ([68]) implies that V is positive 
and it reduces to V — X — 1 for [i ^> A. This simple 
result follows from the fact that recession of the micro- 
tubule is controlled by the unit conversion rate. As soon 
as conversion occurs, detachment occurs immediately for 
fi ^> X and the microtubule recedes by one step. Since 
advancement occurs at rate A, the speed of the tip is sim- 
ply A — 1. However, for extremely large fi the prediction 
V = X — 1 breaks down and the microtubule becomes 
compact. To determine the phase boundary in this limit, 
consider first the case fj, = oo. As shown in the pre- 
vious section, the probability of a catastrophe roughly 
scales as e _7r A / 6 so that the typical time between catas- 
trophes is e" A / 6 . Since V = A — 1, the typical length of 
a microtubule just before a catastrophe is (A — l)e 7r A / 6 . 
Suppose now the detachment rate [i is very large but fi- 
nite. The microtubule is compact if the time to shrink a 
microtubule of length Ae^ A / 6 //x is smaller than the time 
(pA) -1 required to generate a GTP+ by the attachment 
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I ■ • • — ) ==>• I ■ • h) and thereby stop the shrinking. We 

estimate the location of the phase boundary by equating 
the two times to give 

/i, ~ pAV 3 */ 6 when A>1. (70) 

We checked the theoretical predictions and f7D|) 
for the phase boundary numerically (Fig. [5]). For small 
A, the agreement between theory and the simulation is 
excellent. For larger A, the tubule growth is more in- 
termittent and it becomes increasingly difficult to deter- 
mine the phase boundary with precision. Nevertheless, 
the qualitative expectations of our theory remain valid. 

C. Fluctuations of the Tip 

Finally, we study the fluctuations of the tip in the small 
and large A limits. In the former case but also on the 
growth phase pX > (j,, the tip undergoes a biased random 
walk with diffusion coefficient 

D(p,X,fi) = pA + M w hen 1 > pX > n (71) 

For large A, the analysis is simplified by the principle 
that can be summarized by: "The leading behaviors in 
the A — ► oo limit are universal, that is, independent of 
p and /j," This is not true if p is particularly small [like 
A -1 ] and/or if /i is particularly large [like /i* given by 
(|70p ]. But whcnpA 3> 1 and /i -C (J.*, the above principle 
is true, and 

V = X, D=~ (72) 

in the leading order. 

The computation of sub-leading behaviors is more 
challenging. We merely state here two asymptotic re- 
sults. When p <C A, we again have the relation D = V/2, 
with V = A + 1 - p~ l . If /J,* > /i > A > 1, we find 

V = X-1, D=^~ (73) 

The derivation of the latter uses the probabilities X(L,t) 
and Y(L,t) and follows similar steps as in Sect. IIII Bl 

VI. SUMMARY 

We investigated a simple dynamical model of a micro- 
tubule that grows by attachment of a GTP + to its end 
at rate A, irreversible conversion of any GTP + to GDP - 
at rate 1, and detachment of a GDP - from the end of a 
microtubule at rate fi. Remarkably, these simple update 
rules for a one-dimensional system lead to steady growth, 
wild fluctuations, or a steady state. Our model has a 
minimalist formulation and therefore is not meant to ac- 
count for all of the microscopic details of microtubule 



dynamics. Rather, our main goal has been to solve for 
the structural and dynamical properties of this idealized 
microtubule model. Some of the quantities that we de- 
termined, such as island size distributions, have not been 
studied previously. Thus our predictions about the cap 
and island size distributions may help motivate experi- 
mental studies of these features of microtubules. 

A rich phenomenology was found as a function of the 
three fundamental rates in our model. When GTP + at- 
tachment is dominant (A 3> 1) and the attachment is in- 
dependent of the identity of the last monomer on the free 
end of the microtubule (p = 1), the GTP+ and GDP" 
organize into alternating domains, with gradually longer 
GTP + domains and gradually shorter GDP - domains 
toward the tip of the microtubule (Fig. Here, the 
parameter A could naturally be varied experimentally by 
either changing the temperature or the concentration of 
tubulins (free GTP + ) in the solution. 

The basic geometrical features in this growing phase 
of a microtubule can be summarized as: 



symbol 


meaning 


scaling behavior 


N 


# GTP + monomers 


A 


L 


tubule length 


Xt 


(k) 


average cap length 


^/irX/2 


I 


# islands 


A/2 


h 


# GTP+ fc-islands 


2A/fc 3 


Jk 


# GDP - fc-islands 


A/fc 2 



We emphasize that the island size distributions of GTP+ 
and GDP - are robust power laws with respective expo- 
nents of 3 and 2. In the limit of p <C 1, in which attach- 
ment is suppressed when a GDP - is at the free end of the 
microtubule, the average number of GTP + monomers on 
the microtubule asymptotically is pX, while the rest of 
the results in the above table remain robust in the long- 
time limit. 

Conversely, when detachment of GDP - from the end 
of the tubule is dominant (detachment rate fi — > oo, a 
rate that also could be controlled by the temperature), 
the microtubule length remains bounded but its length 
can fluctuate strongly. When the attachment rate is also 
large, the strong competition between attachment and 
detachment leads to wild fluctuations in the microtubule 
length even with steady external conditions. We devel- 
oped a probabilistic approach that shows that the time 
between catastrophes, where the microtubule shrinks to 
zero length, scales exponentially with the attachment 
rate A. Thus a microtubule can grow essentially freely 
for a very long time before undergoing a catastrophe. 

For the more biologically relevant case of intermedi- 
ate parameter values, we extended our theoretical ap- 
proaches to determine basic properties of the tubule, such 
as its rate of growth, fluctuations of the tip around this 
mean growth behavior, and the phase diagram in the 
(A, p) parameter space. In this intermediate regime, nu- 
merical simulations provide more detailed picture of the 
geometrical structure and time history of a microtubule. 
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APPENDIX A: JOINT DISTRIBUTION FOR p = 1 

Introducing the two-variable generating function 

9(x,y,t)= x L y N P(L 1 N 1 t) (Al) 

L>W>0 

we recast (|13[) into a partial differential equation 

^ = ( 1 -f)^-A(l-xy)0' (A2) 

Writing 

9(x, y, t) = e M*y-d-*) in(i-i0] q( Xj y ; t ) (A3) 

we transform (|A2|) into a wave equation for an auxiliary 
function Q(x, y, t) 



8Q , , dQ 

m ={1 - y) dy- 



(A4) 



whose general solution is 

Q(x,y,t)=*(xMl-v)-t) (A5) 

The initial condition P(L,N, 0) = <5l,o<5/v,o implies 
r P{x, y, 0) = 1 and therefore 

e A[*3,-(l-x)ln(l-2/)] $( X) l n (l - y)) = 1 
from which 

$(a,6) = e ^ 1 -°) b +"( eb - 1 )] (A6) 
Combining Eqs. (|A3|) . (|A5|) - (|A6[) we arrive at 

A" 1 In? = xy(l - e _t ) - f - x(l - e _t - t) (A7) 

APPENDIX B: JOINT DISTRIBUTION FOR p 1 

For p ^ 1, we consider the distributions Xn and Yjy, 
defined as the probabilities to have N GTP+ monomers 
with the tip being either GTP + or GDP~, respectively. 
These probabilities satisfy a closed set of coupled equa- 
tions. In the stationary state these equations become 

(N + \)X N = \X N - 1 +p\Y N _ 1 + NX N+1 (Bla) 

(N + p\)Y N =X N+1 + (N + l)Y N+1 . (Bib) 

Since Xq = 0, it is convenient to define the generating 
functions corresponding to Xn and Yj\r a s follows: 



£ z N -'x N 

N>1 



X{z) 
%z) = z N Y N 



(B2a) 
(B2b) 



Now multiply Eq. (|BTa| by z N and Eq. (|B"lb|) by z N ~ 1 
and sum over all N > 1 or N > 0, respectively, to obtain 



P \y = X - C(X - X') 

pAy = x-cy'. 



(B3a) 
(B3b) 



where £ = A(z — 1) and prime denotes a derivative in £. 
We can reduce these two coupled first-order differential 
equations to uncoupled second-order equations: 

CX" + (2+ P A-C)X'-(l+pA)X = (B4a) 
Cy" + (2+pA-C)y'-pAy -0. (B4b) 

The solutions are the confluent hypergeometric functions 

pA 



X{z) 



1+pA 
1 

1+pA 



F(l+pA;2+pA;C) (B5a) 
F(pA;2+pA;C). (B5b) 



These generating functions have seemingly compact 
expressions but one has to keep in mind that the X and 
Y probabilities are actually infinite sums. For instance, 
Yq = y(z = 0) = y(C = — A). Recalling the definition of 
the confluent hypergeometric function we obtain 



Yn = 



1 



1+pA 



F(pA;2+pA;-A) 



1 



E 



(pA)„ (-Ay 



1 + pA ^ (2 + pA)» n\ 



where (a) n = a(a + 1) . . . (a + n — 1) = T(a + n)/T(a) is 
the rising factorial. Note that IIq = Yq. Computing 



1+pA 



F(l+pA;2+pA;-A) 



Y 1 = A 



pA 



(l+pA)(2+pA) 



F(l+pA;3+pA;-A) 



one can determine 11!= Xi + Y\ . Some of these formulas 
can be simplified using the Kummer relation 



F(a;b-X) = e<F(b~a;b;~() 



For instance, 



Y° = E 



n>0 



(n + l)A n e- 
(l+pA)„+ 



x 1 = pA£- 



n „ — A 



A"e 



pA) 



n+l 



(re+ l)A n e 



n„-\ 



Tl>0 



(1 +pA)„ +2 



Af>0 
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